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ABSTRACT* 


The  first  part  of  this  thesis  is  concerned  with  the 
mathematics  of  filtering  in  discrete  time.  Pilters  are  de- 
fined  for  the  purposes  of  1)  condensing  waveforms  into  im¬ 
pulsive  functions*  2)  wave  shaping*  3)  noise  suppression* 

4)  signal  detection  according  to  the  criterion  of  maximum 
signal-to-noise  output  at  an  instant*  and  5)  the  same  over 
an  interval.  The  behavior  of  the  complex  Fourier  trans¬ 
forms  of  some  of  these  filters  is  considered  and  connec¬ 
tion  is  made  with  the  theory  of  orthogonal  polynomials. 

This  leads  to  the  possibility  of  a  feed  back  representa¬ 
tion  of  these  filters. 

In  the  second  part*  computational  experiments  are 
described  in  which  digital  filters  are  applied  to  seismic 
body  waves  to  1)  try  to  determine  whether  the  first  arrl - 
val  is  up  or  down  on  a  seismogram  corrupted  with  micro- 
seismic  noise*  2)  increase  slgnal-to-nolse  ratio  on  seis¬ 
mograms  where  noise  has  almost  obliterated  signal*  3)  as¬ 
sign  polarity  to  each  of  two  seismic  first  motion  wavelets 
so  they  can  be  termed  "same"  or  "opposite*"  4)  remove  spec¬ 
trum  of  seismometer  from  data*  5)  investigate  the  time  var¬ 
ying  spectral  structure  Of  underground  nuclear  shot  seis¬ 
mograms  . 


*  Submitted  to  the  Department  of  Geology  and  Oeophysics  on 
January  14*  196„  in  partial  fulfillment  of  the  require¬ 
ments  for  the  degree  of  Master  of  Science  at  the  Massa¬ 
chusetts  Institute  of  Technology 
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INTRODUCTION 

Although  time  is  a  continuous  parameter  It  often 
happens  that  observations  are  made  at  discrete  time  Inter** 
vals.  Sven  when  continuous  observations  are  made,  it  Is 
often  desirable  to  digitize  them  for  computer  processing. 

This  is  strong  reason  to  do  some  mathematics  In  discrete 
time.  An  even  stronger  reason  as  we  will  see  Is  that 
things  which  are  conceptually  quite  hard  In  continuous 
tiaie  have  analogues  In  discrete  time  which  are  easier  to 
understand. 

Portunately.  In  dlsorete  time  many  general  principles 
can  be  observed  with  wavelets  of  very  short  time  duration. 

This  enables  us  to  consider  some  very  simple  examples  before 
launching  off  into  the  general  theory.  These  simple  examples, 
however,  will  not  forshadow  the  way  in  which  we  will  connect 
the  theory  of  least  squares  filters  with  the  general  theory 
of  orthogonal  polynomials. 

We  denote  time  functions  b^  with  a  subscript  as  the 
time  parameter.  When  the  time  function  has  finite  time 
duration,  we  may  denote  it  as  b  or  (bQ.b1. . . . .bn) .  Any 
time  function  which  has  finite  energy  is  called  a  wavelet. 

The  memory  functions  of  filter*  too.  are  sometimes  called 
wavelets . 


I.  Introductory  Examples 

We  Introduce  the  main  topics  by  means  of  some 
examples.  One  is  given  in  discrete  time  an  Input  series 
(bQ,b1)  of  length  (time  duration)  two.  a  filter  (aQ.a1) 
with  an  impulse  response  of  length  two.  and  the  output 
resulting  from  convolution  to  be  (c0jCi#c2^  °**  length 
necessarily  three.  The  t  Is  determined  from  the*t  and  the 
1>  by  convolution  as  Is  the  usual  procedure  for  linear 
filters,  l.e. 


c 


r  c.  -  a.h* 

C»  *  b# 

.  C*  * 
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or 


1)  Spiking  filter 

To  design  a  spiking  filter  one  would  choose  (a0»ai) 
so  that  3  comes  out  as  closely  as  possible  to  a  spike,  i.e. 
either  (1*0,0)  or  (0,1,0)  or  (0,0,1). 

2)  Wave  shaper 

To  design  a  wave  shaping  filter  one  would  choose 

(a  ,a,)  so  that  3  comes  out  as  closely  as  possible  to  some 
'01  9 

prescribed  waveform  (d0,d^,dg). 

3)  A  matched  filter  « 

To  design  a  matched  filter  one  would  choose  (aQ,a^)  i 

so  that  c^  comes  out  as  large  as  possible  while  making  the 
unit  energy  constraint  ,*s|)  on  the  filter.  In  this 

problem  one  doesn't  care  what  cQ  and  Cg  turn  out  to  be. 

4)  Maximum  energy  sum  filter 

To  design  a  maximum  energy  sum  filter  one  would 
choose  (aQ,a1)  so  that  the  energy  output  (C0  f  &,*"  +  c*  ) 
comes  out  as  large  as  possible  while  making  the  unit  energy 
constraint  (A*'*’  I  )  on  the  filter. 

A  quick  sketch  of  the  solutions  to  these  problems  is 
as  follows:  Since  the  spiking  filter  is  a  special  case  of 
the  wave  shaper  it  will  be  sufficient  to  work  out  the 
solution  for  the  wave  shaper.  Requiring  3  to  be  as  close 
as  possible  to  3  is  equivalent  to  minimizing  the  squared 
distance  between  them 

-J|*  =  (^•"O****'  C*i“4i)  +  (ca“J0  j 

-  («#  i*-  J0 )*V  (*#  ki  +  («,  K 

1 

; 

i 
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Setting  the  partial  derivatives  with  respect  to  aQ  and  a^ 
equal  to  zero  we  get  the  simultaneous  set  for  1^. 


ft +  (b,b#)a, 

(  b,  bp)  6,  +  b,*)tf, 


b9^«  b.  d, 
b9d,  «f  b,  d* 


We  mention  the  particular  case  d« (1,0,0)  called  the 
zero  delay  spiking  filter.  The  solution  of  the  simultaneous 
set  is 

(<*.»<*.)  (k.*+  X*,  * 

Recalling  that  subscripts  are  the  time  variable  we  now 
consider  the  Pourier  transform  of  the  solution 

F*(<*0  =  atf  +  «,e  =  (b9  +  b,  y  -  b„b,e 

The  only  zero  of  this  complex  function  is  in  the  upper 
half  of  the  complex  frequency  plane,  a  fact  which  will  be 
shown  true  for  all  zero  delay  spiking  filters.  This  has 
considerable  importance  in  feedback  systems  and  in  some 
other  connections  to  be  discussed. 

The  solution  to  the  matched  filter  problem  posed  in 
3)  above  is  most  easily  done  by  means  of  Lagrange  multipliers. 
We  wish  to  maximize  c^  under  the  constraint 
Lagrange's  method  is  then  to  maximize 

•*•*■*)] 

Setting  the  derivatives  with  respect  to  aQ  and  a^  equal  to 
zero,  one  gets 
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(«.,<*.) 8  (k->  •>.)/&  - 


Ckjj  i>.L 


Thus  the  filter  (*Qt*±)  is  simply  the  signal  input  time- 
reversed  and  multiplied  by  a  scale  factoi*. 

The  solution  to  the  maximum  energy  sum  problem  4) 
is  somewhat  like  the  matched  filter.  Again  one  uses 
Lagranges  method  and  maximizes 


by  setting  derivatives  with  respect  to  the  components  of  a 
equal  zero.  This  results  in  the  equations 


which  is  the  standard  eigenvector  (a),  eigenvalue  (X) 
problem.  The  two  solutions  to  this  problem  are 


*>w»  «.*  ~  j sr 


It  is  notable  that  the  fourler  transform  of  these  functions 


— 


|f  • 


I- 


-Iu 

e 


have  zeros  on  the  real  frequency  axis.  This  will  also 
happen  with  longer  filters. 

II.  Spiking  Pilters 
A .  Normal  Equations 

In  the  first  Introductory  example  we  considered  the 
problem  of  building  a  two  term  filter  which  would  condense 
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a  two  term  input  into  a  spike  function.  Now  we  would  like 
to  build  an  mtl  term  filter  to  condense  an  n+1  term  input 
into  a  spike. 

•dp 

A  data  wavelet  is  given  by  b«(b0*b1# ...,bn) .  We 
plan  to  construct  a  filter  a  -  (a0,al' • *  *  *aa) •  Filtering 
is  defined  in  this  way:  When  data*S  goes  into  a  filter's, 
an  output  wavelet  t  is  produoed  according  to  the  following 
matrix  multiplication. 


(II-D 


This  operation  is  often  oalled  complete  transient  convolution. 
This  is  more  loosely  written  as 


Here  a  small  amount  of  confusion  can  arise  about  the  limits 
of  the  sunmatlon  because  negative  subscripts  may  appear 
within  the  summation.  What  is  meant  is  that  one  should 
consider  the  terms  "off-the-ends"  of  the  wavelets  to  be  zero. 
With  this  consideration  we  might  write  the  limits  of  the 
summation  as  minus  to  plus  infinity.  The  artifice  of  using 
infinite  limits  on  the  sums  turns  out  to  avoid  some  need¬ 
lessly  cumbersome  notation. 

Now  we  Introduce  another  wavelet  d  which  will  have 
the  same  number  of  components  as  o.  We  call  d  the  desired 
output  of  the  filter.  We  saw  that*£  is  the  actual  output. 
The  actual  output  "t  was  seen  to  be  a  funotlon  of  the  input 
b  and  the  filter  T.  The  problem  now  is  to  determine  ?  so 
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that  "t  and  d  are  very  much  alike.  Specifically  we  will  choose 
"S  80  that  the  difference  vectorfc-S|haB  minimum  length 
squared  (in  nfmfl  dimensional  space).  In  other  words  we 
are  minimizing 


£  CCi  -0* 


(II-3) 


by  varying  the  components  of  3. 
for  “5  in  terms  of  *1  and  Tj  we  get 


iat+fl 


t  *o 


Inserting  the  expression 


(II-*) 


This  function  of  rofl  variables  will  be  minimized  If  its 
partial  derivative  with  respect  to  each  of  (a0*ai» • • . *am) 
equals  zero.  Setting  derivatives  with  respect  to  a^  equal 
zero  we  get  an  expression  for  mfl  equations 


0  =  «.  -<*;) 


(H-5) 


where  one  equation  is  implied  for  each  value  of  fo). 

These  are  called  normal  equations  because  they  say 
that  the  error  vector,  the  quantity  in  brackets,  will  be 
normal  or  perpendicular  to  the  space  spanned  by  the  vector 
set  (c°Iumn  vectors  in  the  matrix  of  equation  II-l) . 

We  bring  the  equations  Into  standard  form  by  bringing 
the  homogeneous  part  (the  part  depending  on  "t)  to  the  left 
side  and  the  inhomogeneous  part  to  the  right 

s  (II-6) 

t*0  fc*0 

In  matrix  form  the  normal  equations  become 


(II-7) 
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which  can  be  abbreviated 


&T(&a)  -  BTd 


(II-8) 


and  whidh  *is  Identically  equal  to 

(BT&)a  -  B  d 

The  matrix  I? B  can  be  written  as 


( II— 9) 


t*  r.  K 
►,  r0  r, 

N 


tm 

h 


where 


t  -0 


0 


if 

if  i 


This  r^  is  called  the  unnormalized  transient 
autocorrelation  of  t>. 
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We  list  three  special  cases  of  these  equations. 

1.  Zero  delay  Inverse  filter  -  This  is  when  *cf  • 

(1,0,0, . • • ,0) . 

2.  Spiking  filter  -  This  is  when  the  impulse  is 
chosen  any  where  in  *3.  It  has  been  frequently  observed 

.  in  practice  that  putting  the  impulse  near  the  middle  of  3 
results  in  an  actual  output  ’S  which  resembles  cl  more 
fclosely  than  if  d  had  been  chosen  as  in  the  zero  delay  case. 

3.  Waveshaping  filter  -  This  is  when  3  is  not  chosen 
to  be  an  impulse  at  all,  but  is  chosen  to  be  some  arbitrary 
wavelet.  The  filter  a  then  tries  to  convert  the  wavelet 

b  into  the  wavelet  *J. 

It  is  worth  noticing  that  the  homogeneous  part  of 
the  normal  equations  (II-9)  depends  only  upon  the  autocor¬ 
relation  of  the  input*?  and  not  on*?  itself.  If  the  desired 
output  of  the  filter  is  an  impulse  with  no  delay  (3  -  (1,0, 
0, .  •  • ,0) )  then  the  inhomogeneous  part  becomes  the  column 
vector 

BtJ  » 


Now  in  this  case  we  see  that  the  waveform  b  does  not  enter 
the  inhomogeneous  part  either,  except  for  the  magnitude  of 
bQ.  Inspection  of  the  normal  equations  shows  that  this 
magnitude  will  not  affect  the  waveform  of  the  filter  "9  except 
as  a  scale  factor. 

Thus  the  normal  equations  in  this  special  case  (zero 
delay  inverse  filter)  depend  upon  the  signal  waveform,  but 
only  through  its  autocorrelation.  Since  autocorrelations 
contain  no  phase  information  it  would  be  a  curious  point  as 
to  what  the  phase  spectrum  will  be  of  the  solution's.  We 
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will  study  this  later  and  come  to  the  curious  conclusion 
that  the  phase  spectrum  Is  such  that  as  much  as  possible  of 
the  energy  In  the  Waveform'S  Is  cramped  up  as  close  as 
possible  to  aQ.  This  Is  called  the  property  of  minimum 
phase  delay  of  the  waveform  a. 

To  fix  ideas  we  now  give  an  example  of  the  deter¬ 
mination  of  a  zero  delay  Inverse  wavelet.  Suppose  that 
the  signal  we  are  dealing  with  Is  the  waveform  b«(2,l). 

We  want  to  design  a  three-term  filter  tU(a0,a1,ag) .  The 
desired  output  must  then  be  n+mfl  •  1+2+1  terms  long  and 
Is  3  ■  (1,0, 0,0).  Prom  (II-10)  rQ  -  5,  r^  -  2,  r2  ■  0. 

The  normal  equations  are 


>  x  6 

T 

3u  r  x 

A, 

- 

0 

0  3-  s_ 

s 

0 

and  the  solution  Is  IT  -  (42,-20,8)/85.  To  see  how  good  the 
filter  Is  we  compare: 

actual  output^  -  (84,2,-4,8)/85 
desired  output  cl  »  (1,0, 0,0) 

B.  Minimum  Phase 

Discussions  of  minimum  phase  In  the  literature  are 
mostly  In  terms  of  continuous  time.  Here  we  wish  to  develop 
its  properties  from  the  point  of  view  of  digital  filters 
which  are  not  so  well  known.  We  begin  by  considering  an 
autocorrelation  function  of  the  type  of  equation  (11-11) 
where 

Y  =  j  *-*+»>  *  (n-13) 

•¥ 

We  wonder  What  functions  b  might  have  this  autocorrelation. 
After  we  have  found  the  class  of  functions  b  that  have  this 
autocorrelation  we  can  enquire  which  one  has  Its  energy  as 
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close  as  possible  to  bQ  and  is,  therefore,  the  minimum  phase 
delay  wavelet.  One  thing  which  we  know  to  begin  with  Is 
that  mors  than  one  wavelet  b  may  have  autocorrelation  r 
(for  example;  the  time  reversed  waveform,  the  negative 
waveform,  and  the  time  reverse  of  the  negative  waveform) . 

We  begin  by  spectral  considerations.  Let  P  denote 
Pourier  transform.  It  is  commonly  known  that  the  energy 
density  spectrum  of  the  wavelet  may  be  expressed  in  two 
equal  ways: 


F>> 


(11-14) 


Thus  the  problem  is  to  factor  Fr(o>)  Into  Fb(fc>)  and  Fb(to). 
Then  we  can  simply  take  the  Inverse  transform  of  Fb(U) 
to  get  the  waveform *£.  The  Pourier  transform  of  r  Is 
simply 


(11-15) 


and  letting  z«  £  we  get 


f=  Ci)  -  'Z'VB+  ■Z'h4\_l+- ••+2  r„ 

~  1  "(hi  ^  ^  (ii- 16) 


•  » 


We  notice  that  the  spectrum  has  been  represented  as  a  poly¬ 
nomial  in  z.  The  usual  procedure  in  factoring  a  polynomial 
is  to  find  its  zeros.  Since  r^-r^,  we  notice  that  Pr(Z)  is 
unchanged  if  we  replace  z  by  z“*.  Thus  if  Pr(ZQ)  is  zero 
then  Pr(l/feQ)  will  also  be  zero.  Thus  for  every  zero  zQ, 
z0“1  is  also  a  zero.  Also  since  the  coefficients  of  the 
polynomial  are  real  the  zeros  are  either  real  or  they  occur 
in  conjugate  pairs.  Thus  if  zQ  is  a  zero  then  ZQ  is  a  zero. 
Most  of  the  zeros  will  probably  occur  then  in  groups  of 
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four  such  as 


15 


where  there  are  two  single  zeros  on  the  unit  circle.  It 
turns  out  that  this  can't  happen.  What  we  are  plotting 
here  is  possible  locations  of  zeros  of  energy  density 
spectra  like  equation  (11-16).  When  z_  ip  on  the  unit 

i  °  tU) 

circle  UJj  is  real  by  the  relation  Z»  &  .  Thus  we 

are  talking  about  the  spectrum  at  some  real  frequency.  A 
function  like  the  following 


which  has  a  single  zero  at  is  not  an  energy  density 
spectrum  because  it  is  not  positive  for  all  frequencies. 

More  generally,  energy  density  spectra  cannot  have  zeros 
of  odd  multiplicity  on  the  unit  circle  in  the .^jplane . 

We  now  know  that  for  every  zerp_fc»  2  C.  0  of  the 

I  I  w 

energy  spectral  polynomial  that  20“*©  *  is  another  zero. 
After  we  factor  the  spectral  polynomial  we  will  be  able 
to  write  the  spectrum  as 

F,  a)  =  z  V'X*-  *;>]  (II_18) 


or  in  terms  of  U) 


A  O) 


6  M  (n-19) 
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Now  If  we  show  A  (<o)  w**  B  (  w)  then  we  have  factored 
the  spectrum  Fr(Uf)  into  the  desired  conjugate  parts 

*b(U>)Fb(U».  ,iu> 

But  both  are  polynomials  In  C  of  order  n  and 

both  A(<**)  and  B(UJ)  have  the  same  zeros.  Thus  they  must 

\ 

be  the  same  function  except  for  a  constant  multiplicative 

factor.  This  can  be  absorbed  from  the  factor  r  c 

n 

This  is  called  factoring  the  spectrum. 

We  notice  that  the  factorization  could  have  been 
done  in  many  ways  depending  on  which  of  the  pair  of  zeros 
is  put  into  Fb(t»0  (the  other  one  then  going  into  «*>)). 
Normally,  there  would  be  2n  different  ways  of  doing  this, 
the  exception  being  the  degenerate  case  when  zeros  occur 
with  multiplicity  greater  than  unity.  Then  there  would 
be  fewer  than  2n  wavelets  with  the  same  energy  density. 

One  of  these  possible  factorizations  is  of  parti¬ 
cular  signif igance .  The  factoring  is  done  so  that  all  of 
the  zeros  which  are  outside*  of  the  unit  circle  are  put 
into  and  the  opposite  member  of  each  pair  which  is 

Inside  the  circle  then  goes  into  Fb  (U?).  In  this  case  the 
wavelet  must  be  real  because  each  root  is  either  real  or 
it  occurs  with  its  complex  conjugate. 

Combining  all  complex  roots  z  with  their  complex 
conjugates^  +t |)we  write  ?or  wavelet's 
transform 

Oh  I  *-*i*  Othfh  IIK«t  out 1 1*4#  CjWc 

Taking  the  Inverse  transforms  and  letting  denote  con¬ 
volution 

1  The  case  with  zeros  exactly  on  the  unit  circle  corresponds 
to  a  spectrum  which  is  exactly  zero  for  some  real  (jJ  .  In 
any  physical  case  one  can  usually  perturb  the  spectrum  slightly 
to  avoid  this  difficulty. 


17 


ta* -  F-*[C  )(  >•••(  X  >— ] 


(n-19.2) 

Thus  we  have  a  string  of  convolutions  of  many  wave¬ 
lets  each  of  either  2  or  3  Instants  duration.  Since  all 
of  the  roo.ts  were  chosen  outside  the  unit  circle  we  have 

|2*|  >1 

This  means  that  in  each  of  the  wavelets  the  first  term  is 
larger  in  absolute  value  than  th>  last.  Thus  in  the 
convolution  of  all  terms,  the  energy  will  be  compacted 
toward  the  beginning.  If  any  one  of  the  zeros  had  been 
chosen  instead,  from  inside  the  circle,  then  the  energy 
would  be  spread  further  out  on  the  time  axis. 

We  will  now  prove  that 

^  in  Int«rv*l  Q  kT 


o£  Wefct 

T  (+/*•) 

the  summed  energy  from  0  to  any  time  t  for  the  minimum 
phase  wavelet  is  greater  or  equal  to  that  of  anyother  wave¬ 
let  with  the  same  spectrum. 
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Consider  a  two  term  wavelet  (b,sy  "bigger,"  "smaller," 
with  its  zero  outside  the  circle.  Convolve  it  into  an 
arbitrary  wavelet  $  »  (po'Pl*  *  *  *  ,P]«)  *  Thc  re8ult  is 

>  *•"♦>.> . 

=  (b  P.>  hr. sii) 

If  instead  we  had  chosen  the  reversed  wavelet  (s,b)  with 
its  zero  inside  the  circle,  we  would  get 

P;n  =  =  + 

Then  we  consider  the  partial  energy  from  time  =  0  up 
to  time  -  T  and  tabulate  the  difference  between  "pin  and  ?out 
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Thus  we  conclude  that  for  any  time  T  the  wavelet  p 

pop 

with  the  zero  outside  the  unit  circle  contains  (b  -s  24 
more  energy  in  the  interval  0*t-T  than  the  wavelet 


out 


with  the  zero  inside.  The  exception  is  at  the  last  lag 


when  they  have  both  put  out  the  same  total  energy.  It  is 


not  difficult  to  show  that  the  above  statements  would  still 


be  true  if  components  of  vectors  were  complex  and  squaring 
were  replaced  with  multiplying  by  conjugates. 

To  prove  the  minimum  phase  wavelet  delays  energy  the 
least,  one  imagines  that  the  convolution  (II-19.2)  had 
been  done  so  that  k  zeros  were  outside  the  circle  and 


n-k  were  inside.  We  have  Just  shown  that  if  one  of  the 
zeros  from  inside  were  replaced  with  an  outside  zero,  that 
the  new  convolution  would  have  less  energy  delay.  This 
argument  is  repeated  until  all  zeros  are  outside. 


Finally,  we  show  that  zero  delay  spiking  wavelets 


determined  by  least  squares  will  have  all  their  zeros  out¬ 


side  the  unit  circle. 

We  recall  the  following  from  previous  portions  of 
this  thesis: 

1)  The  least-squares  spiking  wavelet  is  a  wave¬ 
let  If  which  when  convolved  with  a  given  wavelet  tries 
to  give  an  output  equal  in  the  least  squares  sense  to 
d  ■  (dQ,0,0, . . . ,0) .  Specifically,  is  chosen  to  minimize 


izl 


2)  We  recall  that  the  choice  of  size  of  dQ  affects 
the  solution  vector*^  only  as  a  scale  factor.  Thus  dQ 
could  always  be  chosen  so  that  a0  ■  1.  We  note  that  a 
scale  factor  has  no  effect  on  a  per  cent  total  energy 
graph. 

3)  We  recall  that  if  a  zero  of  a  wavelet  is  removed 
from  inside  the  circle  and  replaced  by  the  conjugate  inverse 
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zero  outside,  that  the  modified  wavelet  has  a  per  cent 
total  energy  curve  which  lies  above  that  of  the  original 
wavelet.  The  per  cent  total  energy  curves  may  touch  one 
another  at  points  except  for  at  time  t-0  where  the  curve 
with  fewest  zeros  inside  the  circle  is  definitely  above. 

We  can  view  the  normal  equations  as  minimizing  the 

dp 

energy  in  a  convolve  b  after  time  t-0  subject  to  the  con- 

p  p 

straint  that  the  energy  at  t-0  be  equal  to  (a0bQ)  -(bQ)  . 

That  is,  we  could  view  the  normal  equations  as  minimizing 
the  per  cent  energy  after  t-0.  But  this  is  the  same  as 
maximizing  the  per  cent  energy  at  t=0.  But  if  the  per¬ 
centage  energy  at  t-0  is  to  be  maximized  for  the  wavelet 
%  convolve  "6,  then  there  must  be  as  few  as  possible  zeros 
inside  the  circle.  This  happens  if *1  has  none  inside  and 
hence  is  minimum  phase. 

C.  Connection  of  Least  Squares  Inverse  Filter  with  Orthogonal 
Polynomials 

Given  an  energy  density  function 

y(a)  =  ... 

one  could  take  that  function  and  use  it  as  a  weighting 
function  to  define  a  set  of  orthogonal  polynomials.  We 
choose  the  interval  of  orthogonality  to  be  the  unit  circle 
in  the  z  plane  which  corresponds  to  the  real  frequency 
axis  from  — VT  to+fr  in  th eh)  plane.  Thus  we  would  con¬ 
struct  a  set  of  polynomials 

-  *,•  +  K  ^ 

In  yL  n 

m  -  n 

(HZ-  ao 


so  that 


s„.  =  {' 
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on  the  real  axis.  Expressing  the  same  thing  with  complex 
polynomials  on  the  unit  circle  one  gets 

s  *£- 

tn*  i  (II-21) 

We  illustrate  the  construction  of  these  polynomials 
in  such  a  way  that  it  will  be  seen  to  be  equivalent  to  the 
least  squares  normal  equations.  Consider  the  construction 
of  fg.  Let  £fn»fm^  denote  the  dot  product  defined  by 
equation  (11-20).  The  vector  fg  is  of  order  two  say 

f*  =  c.  ^  c,i  +  c ^ 

and  must  satisfy  the  orthogonality  conditions 

XM  »® 

ft.  I  I  ]  -0  (11-22) 

‘-I 

Since  fg  perpendicular  to  any  linear  combination  of 
fQ  and  f^>  it  is  perpendicular  to  any  polynomial  of  order 
less  than  2.  Thus  the  orthogonality  relations  could  be 
written 

-  0 

» i  y  *  °  (n-23) 

[f  l1]  1 

This  set  of  orthogonality  requirements  (11-23)  can 
be  written  out  in  full  as 

ij  c.  +  +  £■»  o 

£Mjc.  +  [i.i] c,  -t-LiUfe  o 


(11-24) 
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We  now  examine  the  coefficients  in  this  simultaneous 
set.  Consider  \zn,Z m] 

[l“,  l“]  =  u + JLTjAr*  lui - )  e‘"“« W 

i  r +lr 

(fw  2- U/+ xv-t-'  ■ 


-  K. 


n-m 


Hence  the  orthogonality  relations  (11-24)  can  be 


written 


K  r, 

i-x  h  r.  £, 


CewWrl 


=  O 


(11-25) 


This  is  almost  exactly  the  same  as  the  normal  equa¬ 
tions  (II-7)  for  the  least  squares  inverse  filter.  The 
only  differences  are  a  scale  factor  in  the  inhomogeneous 
part  and  "time  reversal"  of  the  solution.  But  this  will 
not  affect  the  waveform  except  by  a  scale  factor  and 
time  reversal. 

Thus  we  have  shown  the  important  result  that  follow¬ 
ing  two  problems  are  equivalent: 

1)  Find  polynomials  which  are  orthogonal  on  the 
unit  circle  with  weight 

2)  Find  least  squares  zero  delay  spiking  filters 
of  different  lengths  for  the  spectrum 


This  result  is  important  because  it  allows  us  to 
apply  many  results  in  the  classic  field  of  orthogonal 
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polynomials  to  least  squares  filter  theory. 

One  application  is  to  use  the  recurrence  relation  between 
successive  orthogonal  polynomials  to  generate  the  filter  of  length 
n+1  from  the  filter  of  length  n.  This  trick  greatly  facilitates  computing 
the  solution  of  the  normal  equations.  The  relationship  for  getting 
^44  (Z)  from  fm(Z)  is  the  recurrence  relation  (Geronimus  I960) 


&m,o  fm(  V) 


(II-  26) 


where  the  two  side  conditions  used  to  get  *  o  and 


are  first 


f&m.  m  »  (11-27) 


and  second 


(^0*"  “ 


(11-28) 


The  choice  of  sign  for  the  square  root  is  immaterial  as  far  as 

polynomial  orthogonality  is  concerned.  It  is  customary  to  choose  it 

so  thatfl^^  9n+l  first  term  in  the  spiking  wavelet  is  positive.  The 

recurrence  relation  can  be  started  off  by  choosing  any  value  whatever 

The  result  is  just  a  scale  factor  in  the  inverse  wavelet. 

From  equation  (II- 20)  it  is  evident  that  the  recurrence  formula  can  be 

b 

started  off  at  —  CO 

These  relations  appear  to  have  first  been  derived  by  Szego  (1939). 
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Another  readable  account  is  Geronimus  (I960).  Levinson  (1939)  also 
derived  similar,  but  not  identical  relations  for  the  filter  problem*, 
although  he  does  not  mention  any  connection  with  orthogonal  polynomials, 
Levinson's  scheme  is  even  more  useful  than  the  polynomial  recurrence 
relations  because  it  allows  solving  the  normal  equations  for  arbitrary 
inhomogeneous  part. 

Another  valuable  result  of  the  connection  of  filter  and  poly¬ 
nomial  theory  is  the  following.  All  the  zeros  of  all  the  polynomials 
generated  by  the  recurrence  relation  above  are  known  to  lie  inside 


the  unit  circle  (Geronimus  I960).  This  means  that  the  time  reverse  of 
the  associated  filter  is  minimum  phase.  Because  of  this  we  can 


invert  the  wavelet,  i.e.  take  the  inverse  of  its  spectrum. 

— - ; - — =  =  + 


Since  the  polynomial  a(z)  has  no  zeros  inside  the  unit  circle, 


the  infinite  series  b(z)  converges  at  least  up  to  and  including  the  unit 


circle.  This  means  that  the  wavelet  b^  has  finite  energy  and  is  mini¬ 
mum  phase.  The  wavelet  "f  has  a  spectrum  which  is  in  a  least  squares 
sense*  equal  to  1/VO*)*  The  spectrum  of  the  infinitely  long  wavelet"b 
is  exactly  the  inverse  of  the  spectrum  of  *2.  Hence  we  conclude  that 
the  spectrum  oft  is  equal  in  a  least  squares  sense  to  VM  .  Thus 
we  have  found  a  way  to  compute  in  a  least  squares  sense  the  minimum 


phase  wavelet  of  a  given  autocorrelation  function.  Futhermore,  the 
_______________________  w  a 

*  Least  squares  in  the  sense  that  j  ■»  minimized 

where  b(W)  has  power  spectrum  (fits)* 
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computation  is  quite  easy  because  of  the  recurrence  relation.  It  is  the 
most  efficient  method  known  to  the  author  who  has  computed  500  terms 
of  the  minimum  phase  wavelet  in  about  a  minute  on  an  IBM  7090  com¬ 
puting  machine, 

D.  A  Comment  on  Autoregressive  vs  Moving  Average  Representations 
A  question  arises  whether  it  is  more  efficient  to  characterize 
a  stochastic  process  by  the  first  n  terms  in  its  "autoregressive  operator" 
or  by  the  first  n  terms  of  its  "moving  average  operator.  "  What  is  meant 
by  this  is  the  following:  Usually  filtering  is  thought  of  in  terms  of 
convolving  filter  coefficients  la  with  a  data  series.  This  might  be  called 
"moving  weighted  averages"  or  more  commonly,  "moving  averages.  " 
This  is  equivalent  to  multiplying  the  Fourier  transform  of  the  data 
by  that  of  the  filter.  Substituting  z*  •  it  is  equivalent  to  multiply¬ 
ing  z-transform  polynomials  which  convolves  their  coefficients. 

Filtering  could  be  done  in  another  way  called  "autoregression.  "  Instead 
of  multiplying  the  data  polynomial  by  b(z)  one  divides  it  by  the  poly¬ 
nomial  a(z).  This  is  called  "feedback"  filtering  for  reasons  which 
should  be  apparent  to  anyone  who  has  ever  divided  polynomials  by  the 
method  of  synthetic  division  (see  Lanczos  1956). 

By  "efficiency"  we  mean  the  following:  suppose  we  want  a 
filter  to  represent  VWand  it  is  easy  for  us  to  compute  both*£  and  b 
quite  accurately;  in  fact,  we  wish  to  use  many  fewer  terms  than  we  can 
compute.  Which  characterizes  ¥£<•>>  more  accurately  for  small  p, 
the  moving  average  approximation  fc.  f  fc,  *  V  •••  or  the 
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autoregression  approximation  l/ (aQ+a,  z+. *  .  ^-a^z**)  ?  Whittle  (in 

press)  observes  that  the  autoregressive  coefficients  seem  more 
efficient  and  suggests  that  the  reason  is  that  for  the  series  he  deals 
with  (economic),  autoregression  is  a  more  realistic  physical  model* 

The  author  has  also  observed  that  the  autoregressive  coefficients 
seem  more  efficient  in  geophysical  time  series,  but  suggests  a  different 
reason*  When  we  digitize  continuous  functions  we  usually  digitize 
at  a  rate  high  enough  to  avoid  appreciable  frequency  fold  over*  A 
typical  spectrum  looks  like 


these  conditions  apply  a  filter  using  feedback  can  do  a  better  job  for 
the  number  of  components  than  a  filter  which  doesn’t* 
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III  Generalized  Wave  Shaper  with  Noise 

A.  Derivation  of  Normal  Equations 

Here  we  imagine  the  following  model  of  a  physical 

system  to  apply 


Physical  System 


We  want  to  design  a  filter  to  operate  on  the  output  of  the  physical 
system  to  give  us  some  preferred  output.  One  set  of  formulas  will 
enable  us  to  handle  the  following  problems* 

Problem  1*  Given  the  information  wavelet  b^,  the  power  of  the  infor¬ 
mation!  and  the  power  spectrum  of  the  noise,  convert  each  information 
wavelet  b^  which  comes  out  of  the  system  to  some  other  waveform  d^. 
For  example  we  may  be  converting  a  long  drawn  out  function  b^  into  a 
nice  short  one  like  a  spike  or  a  minimum  phase  wavelet.  Of  course,  we 
do  not  want  the  filter  to  respond  very  much  to  noise* 

Problem  2*  Given  the  information  power  spectrum  and  the  noise  power 
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spectrum  design  a  filter  so  that  just  the  information  comes  out  as 
uncorrupted  as  possible*  The  information  might  be  allowed  to  come  out 
with  some  time  delay*  On  the  other  hand  we  might  want  to  predict  the 
information  before  it  comes  out  of  the  physical  system*  To  see  that 
prediction  is  a  reasonable  thing  to  do,  con  side  rJag  an  extreme  case  where 
noise  is  absent9  the  linear  filter  b^  "rings"  for  a  long  while,  and  the 
information  white  light  series  consists  of  impulses  widely  spaced  in 
time.  Of  course  we  cannot  predict  the  onset  of  a  ring,  but  once  a 
ring  starts  we  can  easily  predict  the  rest  of  it* 

Problem  2  was  treated  by  Levinson  and  Problem  1  was  solved 
by  the  author  in  connection  with  some  geophysical  problems.  They 
are  very  little  different.  It  will  be  seen  that  Problem  2  is  a  special 
case  of  Problem  1  so  we  begin  with  Problem  1  and  specialize  the 
results  later. 

Let  b  be  the  signal  wavelet  of  length  n+1. 

Let  A  be  the  optimum  filter  of  length  m+1. 

Let  ct  be  the  desired  convolution  of  a  and  b  of  length  n+m+1. 

Let  l(  be  any  noise  wavelet. 

Let  £  be  a  white  light  series  which  is  convolved  with  u  to  give  a 
statistical  model  of  the  noise  process. 

Let  f*  be  a  white  light  series  of  signal  wavelet  (b)  arrival  times. 

Let  denote  convolution* 

The  input  to  the  filter  is  the  signal  plus  noise,  i.e. ,  (btfytt+U#  {  ). 

The  actual  output  will  be  this  convolved  with  the  filter,  i.e.  (b*/<  #  f  )*a. 
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The  desired  output  is  the  wavelet  d,  occurring  every  time  a  signal 
wavelet  arrives,  i.  e.  t  .  The  expected  sum  squared  error 

is  defined  as:  expected  sum  squared  error  =  expectation  of 


(actual  output- desired  output) 1 


(III-l) 


Since  convolution  is  associative  and  commutative  it  is  valid 
and  will  be  convenient  to  drop  all  asterisks  in  the  expansion  of  the 


above  Square, 


= E  x+  (/*  I A) 


(III- 2) 


By  taking  the  expectation  inside,  it  is  seen  that  the  last  two 
terms  depend  on  E(f  JU  )*  We  will  assume  this  to  be  zero.  This 
means  that  the  signal  wavelets  arrive  at  times  which  are  uncorrelated 
with  the  noise  wavelets. 

We  recollect  the  remaining  terms. 


=  (ba-dJ^E  t/*)**-  («*>*•  Ets*) 


(III- 3) 


From  here  on  the  derivation  will  algebraically  resemble  that  of 
the  spiking  filter.  It  is  convenient  to  rewrite  these  convolutions  in 
subscript  summation  notation#  i.e. 

— ►  1>I 

(b (m‘4) 

(«<0X  — »  a„)  Cu^  a,) 
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Since  we  hooe  to  minimize  the  expected  sum  squared  error  we 
will  take  its  derivative  with  respect  to  each  of  the  independent  variables 
a^  and  set  each  one  equal  to  zero.  Hence 


o 


+  E  U4 .j  Ai 


(HI- 5) 


This  can  be  expressed  in  more  compact  form 

VH  sH-i  =  UA-i  Uk-i 


(III- 6) 


having  noticed  that  R  and  W  thus  defined  are  autocorrelation  matrices 


or  Toeplitz  matrices.  The  expression  simplifies  to: 


0 


~L  *  ,)-£  "  "]  in 


where  £ •  j  is  the  Kronecker  delt  * 

[Utilizing  the  .symmetry  of  the  quantities  in  the  left  hand 
square  brackets  we  can  write; 

[t  La")  K,  +£(5‘)Vi.J]*;=E  l*J  ^  4X 


(III- 8) 
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These  equations  oan  easily  be  rewritten  as  a  matrix 
equation  In  the  same  way  as  with  the  spiking  filters. 

If  the  desired  output  d^  were  Just  the  signal 
possibly  with  some  lag  or  some  negative  lag  (prediction) 
then  the  right  hand  side  no  longer  contains  the  waveform 
but  only  its  autocorrelation.  This  would  be  the 
specialization  to  Levinson's  problem. 

We  now  give  some  examples  writing  equation  (III-8) 
In  matrix  form.^ 

Example  1  The  signal  waveform  b^  -  (2,1).  The  signal 
arrives  with  a  frequency  which  gives  it  an  average  power 
.  The  noise  is  white  and  has  unit  power.  The  filter 
should  have  3  terms.  The  desired  output  is  a  spike  after 
unit  delay,  d  ■  (0,1, 0,0).  The  normal  equations  become 


Example  2  Like  example  1  except  the  desired  output  is 
the  same  as  the  signal  input  with  no  delay.  The  normal 
equations  are  like  example  1  except  the  right  hand  side 
becomes  the  column  vector  (3,2,0)  . 

Example  3  Like  example  2  except  that  the  signal  should 
be  predicted  by  one  time  unit.  The  normal  equations  are 
like  example  1  except  the  right  hand  side  because  the 
column  vector  (2,0,0)*^. 

IV  hatched  Filter 

Suppose  one  is  given  an  autocorrelation  function  of 
a  noise  process  and  also  a,  signal  wavelet.  It  is  desired 
to  detect  the  arrival  of  the  signal  wavelets  in  the  pre-  < 

sence  of  the  noise.  The  method  to  be  used  is  to  filter 
the  incoming  mixture  of  signal  and  noise  and  then  say  that 
signals  arrive  where  there  are  maximums  in  the  output. 
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How  should  the  filter  be  designed?  If  the  noise  were  white 
snd  the  filter  memory  wavelet  had  unit  energy,  then  the 
power  output  of  the  filter  with  noise  as  input  would  be 
unaffected  by  the  frequency  characteristics  of  the  filter* 
Then  the  filter  need  concern  itself  only  with  the  signal. 
Thus  the  Introductory  example  (Section  I,  no.  4)  gives 
the  whole  story  when  the  noise  is  white.  The  result  is 
simply  that  the  signal  filter  coefficients  are  Just  the 
time  reverse  of  the  wavelet  and  the  actual  filtering 
operation  then  amounts  to  crosscorrelation  of  the  signal 
wavdlet  with  the  incoming  data.  If  the  noise  is  not  white 
we  must  do  something  a  bit  more  complicated. 

Using  the  same  notation  as  the  previous  section,  the 
power'  output  of  the  filter  with  noise  input  will  be  the 
quadratic  form  EIS*>  Uj.-/*4e  can  choose  the  magni¬ 
fication  constant  of  the  filter  to  be  such  that  this  power 
is  unity.  This  leads  to  the  constraint 

£  (.S*)  Ug-i  ai  "I  (IV- 

Por  simplicity  we  choose  to  make  the  filter  have 
the  same  length  as  the  signal  wavelet  and  we  choose  to 
have  the  maximum  output  come  when  the  wavelet  is  exactly 
in  the  middle  of  the  filter,  i.e.,  the  nth  lag  of  the 
convolution  where  both  1  and  ?  have  length  n.  Thus  we 
maximize 

(sum  on  l) 

subject  to  the  constraint  equation  (IV-1).  Using  Lagrange 
multipliers  one  maximises 


m<w[  Aitn.i+'XE  Cf  )Uc-* 


(IV-2) 
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We  have  differentiated  terms  exactly  like  this  In  previous 
sections.  Letting *Sp  represent  the  time  reverse  of  the 
signal  wavelet  and  U  represent  the  noise  autocorrelation 
matrix,  we  write  the  result 


o  -  \  +  iXEcnua 


(IV- 3) 


solving  for  1  we  get 


a  =  tr'^XsXEu*) 


(IV-4) 


We  can  usually  Ignore  2  X  B(f  *)  since  It  Just  amounts 
to  a  magnification  factor  In  the  filter. 

In  practice  one  may  prefer  not  to  invert  the  matrix 
In  (lV-4)  or  solve  the  simultaneous  set  (IV-3)  since  there 
Is  an  easy  way  around  It.  One  might  simply  prefilter  the 
data  to  whiten  the  noise  and  then  filter  with  b„.  The 
results  would  be  similar,  the  difference  arising  from  end 
effects. 

More  Is  known  about  the  matched  filter.  Suppose  one 
wants,  to  choose  a  threshold  value  for  the  output  and 
announce  "signal"  whenever  the  threshold  is  exceeded  and 
"no-signal"  when  it  is  not.  Then  one  would  like  to  maximize 
the  probability  of  guessing  correctly.  It  can  be  shown  that 
if  the  noise  is  gaussian,  then  the  matched  filter  and  pro¬ 
per  choice  of  threshold  will  maximise  this  probability. 

V.  Maximum  Energy  Sum  Filter 

Consider  the  following  physical  problem.  A  trans¬ 
ient  signal  waveform  is  sent  through  a  dispersive  media. 

The  media  is  such  that  It  may  badly  disperse  the  wave 
without  altering  Its  spectral  content  a  great  deal.  We 
know  what  spectrum  to  expect  of  the  signal  and  we  know  the 
speotrum  of  the  ambient  noise.  We  would  like  to  design  an 
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apparatus  or  procedure  to  anablo  ua  to  make  a  best  guess 
as  to  when  the  signal  arrives.  The  Matched  filter  is  not 
the  answer  beoause  we  do  not  know  the  exact  signal  shape* 
only  Its  speotrum.  The  spiking  filter  Is  not  appropriate 
for  the  sane  reason.  The  Wlener-Levlnson  filter  tries  to 
make  the  output  look  like  the  signal  Input.  In  this  case 
we  don't  even  know  what  the  input  waveform  should  be*  we 
would  Just  like  to  try  to  decide  approximately  when  It 
arrives . 

A  solution  to  this  problem  Is  to  design  a  filter 
which  puts  out  lots  of  energy  when  the  signal  comes  in 
and  minimum  power  when  only  noise  comes  In.  Thus  our 
decision  would  be  based  on  a  system  like  the  following: 


We  would  search  for  the  time  t  when  the  output  was 

m 

maximum  and  then  we  would  say  that  signal  arrived  between 

time  t_  and  time  t  -T. 
m  n 

Taking  this  model  then*  we  seek  to  maximise 

X  energy  output  of  filter  due  to  signal  in  interval  T 

"  expected  power  output  due  to  noise 

(V-l) 

Notice  the  slsillarity  of  this  problem  to  introductory 
exasg>le  4.  it  will  be  seen  that  it  turns  out  to  be  exactly 
the  same  if  the  noise  is  white. 

Since  we  are  intereated  in  a  computer  application* 
we  again  specialise  ourselves  to  filters  and  algnals  which 
are  disore te  in  tine*  and  speotra  which  whoae  autocorrelations 
are  of  finite  time  duration. 
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Using  the  finite  autocorrelations  of  the  signal  and 
noise  we  define  two  wavelets  b^,  a  signal  wavelet,  and  u^, 
a  noise  wavelet.  This  can  be  done  be  the  procedures  des¬ 
cribed  earlier.  These  two  wavelets  may  have  different 
phase  spectra  than  those  of  our  physical  problem,  but  they 
will  have  the  correct  autocorrelation.  Thus  we  begin  with 
the  definitions  used  earlier: 

a^  m  "ideal”  filter  coefficients  (a^  «  0  if  1<  0  or  i^M) 

b^  -  signal  wavelet  (b^  •  O.lf  i<0  or  i^N) 

u1  -  noise  wavelet  (u^  -  0  if  i<0  or  i>N) 

^ A  -  white  light  series  associated  with  noise  process 
^has  variance  1. 

We  use  subscript  summation  notation;  the  expression 


j  bi-i 


has  an  implied  summation  over  all  values  of  the  repeated 
index  J,  J  goes  from  minus  to  plus  infinity.  Thus  the 
given  expression  is  a  veotor  with  free  index  k  and  IS  the 
complete  transient  convolution  of  a  and  b. 

Expression  (V-l)  for  A  with  this  convention  now 

becomes 


*.l>A-j  aioi-i 


We  notice  that  a  quantity  like  b^_ jb^_^  ia  the 
autocorrelation  matrix  Ba_j  of  the  signal  b^  and  denoting 
likewise  U^_j  -  Uj_i  as  the  autocorrelation  matrix  of  ui# 
the  expression  (V-2)  becomes 


Urn- n 


A 


(V-3) 
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To  try  to  maximize  this  ratio,  we  take  its  partial 
derivatives  with  respect  to  each  of  the  independent  variables 
a ^  and  set  them  equal  to  zero. 


9a 


(V-4) 


Multiplying  by  (Um.n*M*n)  we  get 


(v-5) 


The  derivative  operations  are  the  same  in  each  term, 
working  only  with  the  first  we  get 


9a 


"  (^*  ^4  ^il) 


where is  the  Kronecker  delta.  Now  utilizing  the 
symmetry  of  B^_j  and  the  fact  that  1  and  J  are  dummy  vari- 
ables,  this  becomes 

-  +  &X-i  *i 

-  **-  •,--«*{  (V-6) 

Applying  this  result  in  equation  (V-5),  we  obtain 
0  -  (V-7) 


37 


This  is  the  generalized  eigenvalue  problem.  Further* 
more,  since  B  and  U  are  positive  definite*,  this  problem 

a 

ie  known  to  have  N  diatinot  eigenvector  solutions  for  the 
a1  associated  with  N  eigenvalues  X  m>  The  eigenvalues 
must  be  real  and  positive.  Assuming  that  the  eigenvalues 
are  distinct  we  select  for  our  solution  a^  that  eigenvector 
whloh  is  associated  with  the  maximum  eigenvalue.  We  note 
that  eigenvectors  are  determined  only  to  within  a  scale 
factor.  This  corresponds  to  the  physical  fact  that  the 
energy  power  ratio  (V*l)  will  not  depend  on  the  amplifies* 
tion  of  the  filter. 

Looking  back  to  equation  (V*2),  we  see  that  the  numer* 
ator  is  the  energy  in  the  complete  transient  convolution 
of  a^  and  b^,  and  denominator  is  likewise  for  a^  and  u^. 

The  energy  in  the  convolution  of  two  transients  is  well 
known  to  be  the  integral  of  the  product  of  their  energy 
density  speotra.  Therefore,  if  we  were  able  to  find  • 

i 

another  wavelet  a^  which  had  the  same  amplitude  spectrum 
as  a^,  we  would  have  another  solution  to  our  maximization 
problem. 

Prom  the  z*transform  analysis  described  in  Section 
II,  we  know  that  many  finite  wavelets  may  have  the  same 
spectra.  These  different  wavelets  are  obtained  (by  a 
method  due  to  Wold  and  also  Pejer)  in  the  following  way: 

1)  Compute  the  autocorrelation  of  the  given  wavelet. 

2)  Factor  its  z-transform.  3)  Its  zeros  must  occur  in 
pairs,  specifically  if  Z^  is  a  zero,  then  1/2^  is  a  zero. 

Select  either  one  from  eaoh  pair  and  form  (Z-Z^ (Z-Z2) (Z-Z^) . 

This  is  the  z-transform  of  a  wavelet  with  the  same  auto¬ 
correlation  as  the  given  wavelet.  4)  Normally  there  are 

2n  possible  different  wavelets.  By  the  reasoning  of  the 

preceding  paragraph,  these  should  all  be  solutions  of  our  * 

maximization  problem. 

This  is  an  apparent  contradiction  to  the  fact  that 

the  eigenvalue  problem  (V*7)  it  known  to  poaseaa  a  unique 

Mlo  see  that  B  is  positive  definite  recall  that  ‘Ifli  is  a 
quadratic  fora  representing  the  energy  of  output  when  the 
wavelet  b  goes  into  the  filter  1.  Clearly  this  energy  la 
positive  for  any  real  values  of  a.  This  means  that  B  is 
positive  definite. 
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eigenvector  solution.'  a^  for  the  maximum  eigenvalue  X  max< 

The  contradiction  is  resolved  if  and  only  if  all  of  the  zeros 
of  the  z- transform  of  each  solution  eigenvector  lie  on  the 
unit  oircle.  Then  the  seros  z.  equal  their  inverse  conju- 


and  the  2n  different  selections  of  one  from  each  of  the  n 
pairs  of  zeros  all  generate  the  same  wavelet. 

There  is  a  curious  consequence  of  the  fact  that  the 
zeros  of  the  z -transform  of  this  filter  must  be  on  the  unit 
circle.  It  is  that  the  eigenvectors  must  be  either  symmetric 
or  antisymmetric  (for  example  (2.3,2)  or  (4,0, -4)  respec¬ 
tively).  Whether  it  is  symmetric  or  antisymmetric  depends 
upon  whether  there  are  an  even  or  an  odd  number  of  zeros 
at  the  point  Z«1 . 

This  is  a  simple  consequence  of  the  fact  that  the 
eigenvectors  are  real,  and  any  roots  of  the  z-transform 
which  are  complex  must  occur  in  conjugate  pairs.  By  the 
main  theorem,  they  must  also  lie  on  the  unit  circle.  Por 
the  root  s  Cd  «*•  i  we  may  then  state 

«*■+  i 

and 

(a »(2  -*-Lp)(a-«(+ 

=  £*■-  a_cti  +(«*'»■  <»*) 

-  2X  -  1**  +■  I 


Hence,  the  coefficients  of  the  second  order  and  the  zero 
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order  terns  In  z  are  identical  for  all  and  (3  and  the 
wavelet  is  symmetric .  The  same  is  evidently  true  for  all 
the  oomplex  roots.  The  net  convolution  of  all  these 
symmetric  wavelets  is  symmetric.  Hence,  the  eigenvector 
would  have  to  be  symmetric  if  all  the  zeros  were  complex. 
However,  we  also  have  the  possibility  of  zeros  at  two 
places  on  the  real  axis,  *1,  and  +1 .  The  -1  corresponds 
to  symmetric  wavelet  (1,1),  and  the  +1  corresponds  to  the 
antisymmetric  wavelet  (-1,1).  Convolution  by  the  first 
leaves  the  eigenvector  symmetric,  but  an  odd  number  of  con¬ 
volutions  by  the  second  leave  the  eigenvector  antisymmetric . 


Numerical  Example 


Let 

bi.  C  (»,-■,  a.) 
«i  ■  *  (  %  a. ,  i ) 

and 


\\ 

10 

*f 

1h 

-s  C 

u  = 

10 

to 

B  = 

-s* 

/H  -S’ 

to 

XI 

J 

-S’  IH 

solving 

[b  -  xu)«. « •> 


we  get 

Xt  =a.7J.H3L 


X 

X 


X. 

> 


'a  -  ( 

-0 

’a  =  U  j  0 


<n  ot 

.am 
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The  eigen-values  are  distinct.  The  eigenvector 
solutions  for  the  maximum  and  minimum  eigen-values  are  seen 
to  be  symmetric,  and  the  remaining  eigen- value  has  an  anti¬ 
symmetric  eigen-vector.  The  zeros  of  the  z- transform  of 
the  elguh-vectors  are  then  computed  and  plotted: 

A 

The  magnitudes  of  all  the  zeros  are  seen  to  be 
equal  to  1 


B.  Maximum  Energy 'Sum  Pilter  from  Spectral  Considerations 
We  consider  the  same  problem  of  determining  a  filter 
of  finite  length  in  discreet  time  which  is  optimum  in 
the  sense  that  it  maximises  the  ratio: 
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(energy  output  of  filter  due  to  signal) 
(expected  power  output  due  to  noise) 


This  tine  we  solve  it  in  the  frequency  domain  rather 
than  the  time  domain.  Define  the  filter  energy  spectrum 
as  A (to),  the  signal  energy  spectrum  as  B(fc>)  and  the 
noise  power  speotrum  as  0(U>).  Then  the  above  ratio  may 
be  written: 


5*  AO)G<W*« 


If  the  maximum  of  this  ratio  is  finite  then  it  is 
necessary  that  for  perturbations  in  A (to)  we  will  have  I  -  0. 

Since  A(CJ )  appears  in  both  numerator  and  denominator 
it  is  clear  that  a  multiplicative  scale  factor  in  A(&>) 
will  be  unimportant,  in  other  words  we  can  choose  the  scale 
factor  as  we  wish.  In  fact*  we  can  choose  it  so  that  the 
Integral  in  the  demon lna tor  is  some  constant,  i.e. 


|  =  J’*/"  AC«)GOo)A«* 


(V.B.2) 


Then  the  problem  can  be  restated  as  maximizing  the 
numerator 


£>  =  ACM  BMfi, 


(V.B.3) 


subject  to  the  constraint  equation  (2) . 

This  is  a  classic  problem  in  the  calculus  of  varia¬ 
tions  (see  for  example  Hildebrand,  Methods  of  Applied  Math, 
section  2.6).  The  procedure  is  to  maximise  the  quantity 


j '  a  (  At*)  -  \j  At*) tt*)  A* 


(V.B.4) 


subject  to  no  oonstraint.  And  than  later  X  oan  be 
determined  by  (2).  X  i®  oalled  the  Lagrange  multiplier. 
Thus  we  eolve  the  problem: 


O  a  fjj^C Ot<»> -*&<*))  A  «•>*<• 


(V.B.5) 


Since  we  are  dealing  with  functions  In  discrete  time# 
the  spectra  in  equation  (5)  will  all  be  periodic  with 
period  2  (Nyquists).  The  spectra  are  also  even  functions 
of  Cil  .  Therefore#  A#  B#  and  0  oan  always  be  written  as 


AO»)  =  ct,  +  \  ■*<•> 

HM 

so 

-  f^e  <*•»  0-*.nu> 

ml 

GM  ~  ^o+3“X  C++.r\<* 

e»l 


(V.B.6) 


Pourier  cosine  series  whose  coefficients#  the  Greek  letters# 
can  be  recognized  as  the  autocorrelation  functions  of  the 
respective  time  functions.  The  limit  on  the  summation  for 


A(U>)  is  finite  because  the  filter  a^  was  choaen  to  have 
finite  length  and  hence  so  must  its  autocorrelation.  We 
apply  these  forms  to  equation  (5). 
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The  variation  ia  intended  to  be  over  the  correlation 
function  of  the  filter  impulse  response  a^.  The  OC £ 

are  not,  however,  allowed  to  be  varied  arbitrarily,  they 
must  only  be  varied  in  such  a  way  as  to  keep  the  energy 

density  h(OJ)  positive  for  all  Id  .  In  other  words  an 

arbitrary  seleotion  of  the  numbers  0t  i  may  not  really  be 
an  autocorrelation  function.  Therefore,  we  will  express 
the  0(  i  in  terms  of  the  Impulse  response  a^  and  do  the 
variation  in  tenu  of  the  a^  instead,  because  any  set  of 
numbers  a1  is  a  valid  impulse  response.  The  expressions 
relating  ot^  end  a^  are: 

«0  -  «o +  •! +  4 + . ♦  4 

&  1  ■  a0al  +  ala2  +  a2a3  •  *  *  *  +  aH-laN 

to  2  "  a0a2  +  ala3  +  a2a4  +  •  •  •  +  aH-2aH 


to N "  a0aH 


(V.B.8) 


Performing  the  variation  merely  amounts  to  writing 
the  Euler  equations  in  terms  of  the  a^,  the  being 
complstely  independent  variables.  Our  integral  is  of  a 
particularly  simple  form,  therefore,  we  can  obtain  greater 
insight  by  performing  the  integration  directly.  Then  we 
can  set  the  variations  (derivatives)  with  respect  to  the 
a^  equal  to  sero. 

The  integrand  is  the  produot  of  two  cosine  series. 
Using  the  orthononfiallty  of  these  cosines  over  the  interval 
♦  Yr  to  -  fr  equation  (7)  becomes  on  integration 

0  2  f(Pe“V*i<Jcc0‘f**.5[  (P*» 


(V.B.9) 
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* 


It  is  noteable  that  the  formulas  (9)  no  longer 
contain  the  Infinite  sum  which  is  in  formula  (7) .  This 
Important  result  will  be  referred  to  later,  it  means 
that  only  N  lags  of  the  signal  and  noise  autocorrelations 
are  needed  for  the  solution,  Mfl  being  the  length  of  the 
impulse  response  of  the  filter  a^  which  we  are  construct¬ 
ing. 

We  now  differentiate  the  OC^  in  equation  (8)  with 
respect  to  the  independent  variables  a^.  This  may  be 
wri tten: 


*J-n  +  aJ+n 

where  o£j  £n 
Oin  -N 


and  a^O  if  i  0 

aiS.°  if  i  H 

We  now  insert  this  into  formula  (9)  and  reorder 
terms  according  to  increasing  subscripts  of  aj.  This 
step,  although  it  is  complicated  amounts  to  straight¬ 
forward  symbol  manipulation.  The  final  result  can  be 
written  as  the  following  matrix  equation: 


Thus  we  are  led  to  the  same  result  as  the  time  domain 
considerations  (Y-7)»  One  wonders  whether  there  might 
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be  a  useful  connection  here  with  the  general  theory  of 
eigenfunctions  as  there  were  useful  results  of  connecting 
least  squares  filters  with  the  theory  of  orthogonal 
polynomials . 

It  Is  possible  and  seems  likely  that  some  of  the 
statements  about  decision  rules,  maximum  likelihood,  etc. 
which  are  made  about  matched  filters  in  Gaussian  noise 
could  also  apply  to  the  maximum  energy  sum  filter*. 

This  Is  a  topic  which  does  not  appear  to  have  been  Inves¬ 
tigated. 


*  This  possibility  was  suggested  to  the  author  by  both 
Professor  E.  M.  Hof s tetter  and  Professor  T.  R.  Madden. 
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SECTION  VI  First  Motion  Spiking 

A.  Object  and  Motivation 

The  direction  of  first  motion  of  the  ground  at  a 
seismic  station  has  received  considerable  attention  in 
nuclear  detection.  The  essential  idea  is  that  the  first 
motion  resulting  from  an  explosive  blast  should  always 
be  upward  and  away  from  the  epicenter  while  this  would 
probably  not  be  true  for  more  than  half  of  the  time  for 
naturally  occurring  seismic  events.  This  criterion  has 
been  shown  to  be  a  reasonable  one  for  the  Logan  and 
Blanca  test  shots  for  distances  less  than  about  700  km 
(Romney,  1959) .  The  primary  difficulty  in  considering 
seismograms  taken  at  greater  distances  was  the  reduced 
slgnal-to-noise  ratio  further  aggravated  by  the  fact 
that  the  first  motion  was  in  almost  all  cases  smaller 
than  the  immediately  following  oscillations.  On  some  of 
the  seismograms  taken  at  greater  distances  the  first 
motion  appeared  to  be  in  the  wrong  direction  despite  a 
fairly  strong  signal-to-noise  ratio.  The  motivation  of 
the  experiment  to  be  discussed  is  that  perhaps  the  oscilla¬ 
tions  immediately  following  the  first  motion  also  contain 
information  about  the  polarity  of  the  first  motion,  but 
contain  this  information  in  some  latent  way.  This  idea 
is  not  new,  but  no  effective  method  has  yet  been  applied 
to  extract  this  information. 

A  mathematical  technique  for  extracting  this  type 
of  information  is  the  spiking  filter. 

B.  Method  and  Philosophy  of  the  Experiment 

First  a  wavelet,  the  first  motion  and  several  sub¬ 
sequent  wiggles,  is  selected  from  a  relatively  near-shot, 
low-noise,  seismogram.  Then  a  filter  is  designed  such 
that  with  the  wavelet  as  input,  it  will  produce  little 
or  no  output  before  and  while  the  wavelet  is  entering  the 
filter,  a  large  positive  spike  when  the  wavelet  has  fully 
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entered  the  filter,  and  little  or  no  output  thereafter. 
The  filter  la  also  designed  to  have  little  output  when 
naturally  occurring  microselsms  are  its  only  Input.  In 
practically  all  cases,  a  filter  cannot  be  designed  to  do 
these  simultaneous  tasks  exactly,  but  the  one  designed 
does  them  in  the  least-squared-error  sense.  That  the 
ultimate  error  will  be  sufficiently  small  for  practical 
purposes  must  be  tested  computationally. 

The  filter  is  then  applied  to  a  seismogram  with 
a  poorer  signal-to-noise  ratio  which  may  be  at  a  differ¬ 
ent  orientation  to  the  seismic  event  and  at  a  greater 
distance.  If  the  filtered  seismogram  consists  of  low 
levels  noise  preceding  the  abrupt  arrival  of  a  spike  of 
positive  polarity  we  might  then  infer  that  the  direction 
of  first  motion  is  the  same  at  the  second  station  as  it 
was  at  the  first.  If  the  impulse  had  negative  polarity 
we  would  infer  that  the  second  signal  had  undergone  a 
180*  phase  shift  with  respect  to  the  first  signal.  If 
no  impulse  showed  clearly  through  the  background  noise, 
we  would  infer  that  this  experiment  waB  not  successful. 

To  be  more  precise,  in  least-squares  fitting  to 
a  positive  impulse  we  are  assigning  a  polarity  to  a 
clear  first  arrival  wavelet;  then  we  produce  a  filter 
which  can  be  applied  to  wavelets  from  other  seismograms 
of  the  same  event  which  assigns  a  polarity  to  each  of 
these . 

Finally,  we  are  in  a  position  to  examine  the 
possibility  that  the  polarity  is  the  same  at  all  orien¬ 
tations  from  the  source.  If  it  is,  we  infer  that  the 
source  has  rotational  symmetry  and  is  probably  not  of 
natural  origin.  If  the  polarity  on  the  first  clear 
arrival  wavelet  is  assigned  according  to  the  direction  of 
first  motion,  and  if  wiggles  subsequent  to  the  first 
motion  really  do  contain  latent  information  about  the 
first  motion,  then  the  hypothesis  tested  by  this  experi¬ 
ment  is  very  similar  to,  although  not  exactly  the  same 
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as,  the  hypothesis  that  the  first  motion  caused  by  a 
nuclear  explosion  must  be  up  and  away  at  all  source  orien¬ 
tations.  To  point  out  this  difference  more  clearly,  con¬ 
sider  the  seismograms  mentioned  earlier  on  which  the  first 
motions  appeared  to  be  in  the  wrong  direction.  Possibly 
the  first  motion  was  in  the  right  direction  and  obscured 
by  the  noise,  but  it  might  actually  have  been  in  the 
wrong  dlreotion.  Even  if  it  was,  its  polarity  as  deter¬ 
mined  by  the  first  few  wiggles  might  have  been  the  same 
as  that  of  other  seismograms  of  the  same  nuclear  event. 

C.  Choice  of  Parameters 

Several  of  the  seismic  records  from  the  Logan  under¬ 
ground  nuclear  explosion  were  picked  by  eye,  that  is,  the 
first  motions  were  identified  approximately  and  the  first 
3.3  to  4.0  seconds  of  the  seismic  trace  were  considered 
to  1?e  ithe  essence  of  the  signal  wavelet.  The  section  was 
then  tapered  smoothly  to  zero  on  each  end.  The  exact 
way  in  which  this  was  done  ia  depicted  in  Figure  1.  Only 
the  shorter  of  the  two  wavelets  shown  (the  bottom  in  each 
frame)  was  used.  The  wavelet  length,  about  3*73  seconds, 
was  selected  because  it  is  long  enough  to  include  the 
requisite  "first  few  wiggles"  but  not  so  long  as  to  make 
the  solution  of  the  simultaneous  equations  excessively 
time  consuming.  A  sixty  point  inverse  wavelet  which  is 
three  seconds  in  length  at  our  standard  digitization  rate 
requires  about  one  minute  of  IBM  709  time  to  compute. 

The  choice  of  a  method  of  tapering  the  ends  of  the 
wavelet  was  rather  arbitrary.  It  was  motivated  by  two 
eonsiderationat  1)  The  time  of  the  first  motion  arrival 
could  not  be  determined  exactly,  and  to  be  sure  that  the 
first  motion  arrival  was  lnoluded,  about  3/U  second  of  the 
seismic  trace  before  the  apparent  arrival  was  included  in 
the  wavalet.  Sinoe  It  was  also  felt  that  the  wiggles 
nearest  the  first  motion  probably  contain  the  isost  infor¬ 
mation,  wiggles  further  away  were  also  tapered  in  amplitude. 


49 


2)  If  the  wavelet  were  Just  extracted  from  the  seismogram. 

It  would  be  likely  that  there  would  be  strong  discontin¬ 
uities  in  both  the  function  and  its  derivatives  at  these 
ends.  It  would  be  undesirabl*  if  the  spiking  filter 
turned  out  to  be  particularly  sensitive  to  these  artifi¬ 
cially  caused  discontinuities;  hence,  they,  too,  were 
removed  by  tapering. 

To  select  the  coefficients  of  the  spiking  filter, 
the  following  quantity  was  minimized: 

sum  of  square  error  * 

(delta  function  minus  the  convolution  of  the  filter 
with  the  wavelet)2 

O 

♦  2 (the  convolution  of  the  filter  and  the  noise) 

The  noise  referred  to  in  this  expression  is  the 
mlcroaelsmic  noise  which  Just  preceded  the  arrival  of  the 
signal  wavelet.  The  2  in  the  second  term  on  the  right  in 
the  above  expression  was  selected  on  the  basis  of  results 
of  earlier  crude  computational  experiments.  The  choice 
of  the  delay  in  the  delta  function  in  the  first  term  on 
the  right  in  the  above  expression  was  made  such  that  the 
filter  would  be  acting  on  all  of  the  terms  in  the  wavelet 
at  the  time  of  the  filter's  spike  output. 

The  length  of  the  spiking  filter  was  chosen  to  be 
equal  to  the  length  of  the  wavelet,  not  because  of  theoret¬ 
ical  necessity,  but  because  it  was  thought,  for  various 
reasons,  to  be  a  reasonable  choice. 

The  choice  of  practically  all  of  the  parameters  in 
the  above  discussion  is  somewhat  arbitrary.  They  were  all 
selected  initially  on  Intuitive  grounds.  Some  have  been 
more  or  less  Justified  by  almnle  couputational  teats,  others 
remain  to  be  investigated. 

D.  Results 

Aa  a  oheck  on  the  computations  and  a  check  that  the 
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sum-square-error  would  be  small  enough  to  make  the  scheme 
useful,  the  spiking  filter  was  applied  to  the  seismogram 
from  which  it  was  derived.  This  is  presented  in  the  upper 
left  and  lower  right  frames  on  Figure  2  and  Figure  3.  It 
is  seen  that  the  noise  preceding  the  first  motion  is  in 
all  cases  reduced  and  that  the  first  motion  is  condensed 
to  a  neat  spike.  Just  as  it  should  be.  The  "hash"  which 
is  near  the  beginning  and  end  of  some  of  the.  convolution 
traces  is  the  result  of  applying  a  filter  onto  the  ends  of 
a  finite  segment  of  data. 

The  conclusion  to  be  drawn  from  the  first  part  of 
the  experiment  is  that  a  least  squares  error  filter  can  be 
determined  with  the  resulting  error  small  enough  that  it 
will  be  useful  in  simultaneously  reducing  noise  energy 
and  condensing  a  particular  waveform  into  a  spike. 

The  next  part  of  the  experiment  was  to  apply  these 
spiking  filters  to  other  seismograms.  The  spikes  still 
seem  to  be  present  although  they  are  almost  down  to  the 
level  of  the  noise.  This  is  shown  in  the  lower  left  and 
upper  right  frames  in  Figure  2  and  Figure  3.  In  some 
cases  the  noise  before  the  first  motion  appears  to  have 
increased  after  filtering.  This  is  because  all  of  the 
traces  on  the  figures  were  scaled  to  have  a  certain  maxi¬ 
mum  amplitude  amenable  to  scope  display.  Since  the  spike 
was  always  smaller  in  cases  when  the  spiking  filter  was 
applied  to  other  records,  the  resulting  displays  were 
amplified.  The  spikes  generated  from  the  application 
of  spiking  filters  upon  other  seismograms  are  not  clearly 
distinguishable  from  the  noise  in  all  cases.  The  conclu¬ 
sion  to  be  drawn  from  this  is  that  the  first  motion 
wavelet  loses  much,  but  not  all,  of  its  character  in  going 
from  the  station  at  l8G0  km  to  the  station  at  1900  km. 

This  must  be  qualified,  however,  for  the  loss  of  character 
might  not  be  quite  as  great  as  it  would  first  appear;  it 
should  be  remembered  that  the  wavelet  as  determined  at  one 
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station  also  includes  the  noise  at  that  station,  hence  even 
if  there  were  no  change  in  the  wavelet  at  all  during  trans¬ 
mission  from  one  station  to  the  next,  there  will  be  a 
double  corrupting  effect  in  this  computation  due  to  the 
different  noise  at  the  two  stations  which  cannot  be  com¬ 
pletely  eliminated. 

7.  Possible  Modification  to  and  Experimentation  on  this 
Mathematical  Technique 

The  operation  of  the  spiking  filter  in  this  experi¬ 
ment  had  the  undesirable  effect  of  increasing  the  high  fre¬ 
quency  noise.  As  a  result  of  this,  the  filtered  data  looks 
much  more  spiky  than  the  unfiltered  data  making  it  more 
difficult  to  observe  a  true  spike  in  the  filtered  data. 
Heuristically  the  reasons  for  this  are  as  follows.  The 
energy  in  the  spectra  of  the  signal  wavelet  and  the  noise 
tends  to  be  primarily  at  low  frequencies.  If  we  were 
ignoring  noise  and  considering  an  infinitely  long  Inverse 
wavelet,  its  spectrum  would  be  Just  the  inverse  of  the 
spectrum  of  the  signal  wavelet  and  in  this  case  would  con¬ 
tain  very  high  frequencies.  Since  the  filter  is  also  expected 
to  reject  noise  of  low  frequency,  the  result  is  a  filter 
which  is  very  sensitive  to  high  frequencies  and  hence  high 
frequency  noise.  An  important  conclusion  of  this  experiment 
is  that  something  should  also  be  done  about  high  frequency 
noise.  The  analysis  suggests  how  to  make  the  filter  insen¬ 
sitive  to  any  type  of  noise  of  known  autocorrelation. 

Another  approach  is  not  to  require  an  impulse  to  be  the 
output  of  the  filter,  but  instead,  some  wider  burst.  Reason¬ 
ing  again  from  the  limiting  case  of  filters  and  signals  of 
infinite  extent,  this  would  be  advantageous  because  the  pro¬ 
duct  of  the  spectra  of  the  filter  with  that  of  the  wavelet 
must  equal  the  spectra  of  the  desired  response.  B|y  desiring 
a  response  of  a  wide  burst  Instead  of  a  spike  we  may  expect 
to  get  a  filter  less  sensitive  to  high  frequency  noise. 
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0.  Nor*  Possible  Applies tions  to  Nuclear  Detection 

One  could  try  the  following  different,  though  similar 
experiments! 

1)  On  records  taken  at  the  same  distance  and  at  the 
same  station  try  filters  generated  from  a  wavelet  from  one 
nuclear  event  on  a  seismogram  from  another  nuclear  event. 

2)  On  a  record  with  a  clear  first  motion,  compute 
the  spiking  filter  and  then  oonvolve  the  whole  reoord  with 
it,  in  searoh  for  later  arrivals  of  the  same  waveform.  (If 
later  arrivals  are  detected  their  time  delays  can  be  deter* 
mined  to  the  accuracy  of  the  digitalisation  sampling.  Slnoe 
this  is  1/feO  of  a  second,  it  may  lead  to  Improvements  in 
depth  determination  aoouraoies.) 

SECTION  VII  Prediction  Error  Experiment 
I .  Philosophy 

Mioroseismic  noise  can  be  predicted.  Por  example, 
it  was  found  that  given  past  values  on  our  noise  seismograms, 
one  can  easily  predict  1/4.0  of  a  second  into  the  future 
with  an  error  in  power  of  less  than  5#.  Suppose  we  form  a 
new  signal  by  subtracting  the  predicted  seismogram  from  the 
actual  seismogram.  This  new  signal  is  called  the  prediction 
error  signal.  The  amplitude  of  the  prediction  error  signal 
is  expected  to  be  small.  If,  however,  at  some  time  in  the 
mioroseismic  trace  a  real  signal  arrives,  it  cannot,  of 
course,  be  predicted  from  the  noise.  Hence,  at  that  tim 
the  prediction  error  signal  should  suddenly  attain  a  large 
amplitude.  Por  example,  during  the  digitisation  of  our 
seismogram  one  of  the  timing  marks  was  accidentally  traced. 
Naturally  the  timing  mark  oould  not  be  predicted  on  the 
basis  of  the  nolae  whloh  preceded  it.  The  result  was  a 
large  prediction  error  at  that  time.  This  is  deploted  in 
Pigure  1. 

The  mathematical  theory  of  predicting  stationary  time 
series  at  unit  prediction  distanoe  also  showa  that  the  pre¬ 
diction  error  of  a  pure  noise  signal  will  be  a  white-light 
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series.  The  arrival  of  a  signal,  if  it  has  a  different 
speotrum  than  the  noise,  will  result  in  non-white  series. 

Thus  a  person  attempting  to  find  a  seismic  signal 
arrival  by  examining  the  prediction  error  will  look  for: 

1)  large  inorease  in  amplitude 

2)  change  in  white  character  of  trace. 

There  is  another  peculiarity  of  the  prediction  error 
signal.  Its  power  spectrum  is  Independent  of  the  seismometer 
and  recording  system.  This  is  true  both  before  and  after 
P  wave  arrival.  Before,  the  spectrum  is  simply  white.  After, 
it  is  a  function  only  of  earth  motion  power  Bpectra. 

Another  property  of  the  prediction  error  trace  is  that 
the  ratio  of  power  after  to  power  before  signal  arrival  must 
be  anrimprovement  from  the  original  seismogram. 

All  of  these  properties  w'.ll  now  be  derived. 

II.  Mathematical  Derivation 

The  concept  of  prediction  is  treated  in  greater 
detail  elsewhere  (Robinson  195*0  •  The  formulas  are  briefly 
derived  here  in  an  heuristic  manner. 

First  we  make  the  following  definitions.  Let 
s  be  the  given  stationary  series 
w  be  the  one  sided  wavelet  with  the  same 
spectrum  as  s.  of  length  n 
x  be  the  white  light  series  which  when 
convolved  with  w  gives  s 
v  be  the  wavelet  which  is  inverse  to  w  of 
length  m 

d  be  the  predicted  s  at  d  time  intervals 
in  the  future 

m  or  n  or  both  may  be  infinite. 

Let  negative  aubaorlpts  refer  to  the  past,  the  zero  sub¬ 
script  to  the  present,  and  positive  subscripts  refer  to 
the  future.  Let  denote  convolution.  The  white  light 
series  x  can  be  generated  for  all  past  time,  up  to  and 
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Including  the  present  instant  by  the  convolution  or  s  with 
vj  i.e.. 


* 


-  s*v 


(1) 


The  white  light  series  corresponds  to  arrival  times 
of  the  wavelet  w.  The  situation  is  depicted  in  the  sketch 
below. 


To  find  the  predicted  value  of  the  series  at  the 
time  (nowfd)  we  sum  up  the  effect  of  all  wavelets  arriving 
in  the  past.  Those  wavelets  which  may  arrive  between  now 
and  the  time  we  are  predicting  will  contribute  to  the 
error  of  the  prediction. 

Referring  to  the  sketch  above,  our  prediction  now 
pQ,  for  the  value  of  s^  at  the  future  time  l«d  is  thus 
written: 


Fm  (f)  s  XftMk  ^ ^  t 


1*0 


More  generally,  the  prediction  Pj(d)  for  the  value 
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ol  the  eeriee  Bj+(J  at  time  J+d  le  written: 

We  write  this  symbolically  as 

Pfi)  » 

where  w^  ie  the  wavelet  w,  truncated  of  its  first  d  terms* 
Utilising  (1)  and  the  commutlvity  of  convolutions  we  get 

POO  a  S*(v*v4r) 

and  we  can  Identify  v*wt  as  the  predicting  filter. 

The  prediction  error  e^  is  defined  as 

ej(d)  a  actual  series  -  prediction  of  series. 

It  is  a  function  of  the  prediction  distance  d.  We 
will  now  show  that  if  dal*  the  operator  which  generates 
Sj(l)  from  Bj  takes  on  a  particularly  simple  form,  and  ej(l) 
must  be  a  white  light  series. 

Denoting  z- transforms  by  capital  letters*  the  z- trans¬ 
form  of  the  truncation  of  wt  corresponding  to  d-1  is 

¥/(*)-  W# 

The  z- transform  of  the  predicting  filter  is  then 

The  s-transform  of  the  prediction  error  filter  is 
Just  - 

l  -vWIvjW-Hi) 

a  \  -  vet)  w(^+v«vC0 
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But  w  and  v  are  inverses  and  also,  wQ  -  l/vQ,  henoe 
the  s- transform  of  the  prediction  error  filter  is 

*1  -1+  x  v<A 

-■  >A  v(r) 

Hence,  the  prediction  error  filter  is  Just  the 
inverse  wavelet,  scaled  so  that  the  first  term  of  the 
filter  is  +1.  Since  the  prediction  error  filter  is  the 
inverse  to  the  wavelet  of  the  stationary  series,  it  must 
whiten  the  series. 

What  happens  to  the  spectrum  if  a  signal  arrives 
somewhere  on  the  noise  record?  Letting  S  denote  spectrum, 
the  condition  that  the  noise  be  whitened  is: 

S  (earth  noise)  S  (seism,  system)  S  (prediction  error  filter)  »  1 

The  spectrum  of  our  final  graph  is  then: 

S (graph)  -  S  (earth  signal  +  noise)  S  (seism,  system)  S  (pre¬ 
diction  error  filter) 

Combining  the  above  two  expressions  we  get: 

„  / _ S  (earth  nignal  +  noise) 

S  (graph)  -  fc  " (earth  *'T5xse) - 

which  is  independent  of  the  transfer  function  of  the  seismo¬ 
graph. 

This  derivation  contains  some  hidden  mathematical 
aaSumptlons  which  should  be  valid  in  any  real  case . 

(Seismograph  system  is  linear  and  dissipative.  Ground  motion 
satisfies  Paley-Wiener  criterion.) 

The  proof  that  the  prediction  error  filter  must 
improve  the  signal-to-noise  ratio  is  omitted.  It  is  based 
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on  the  fact  that  the  prediction  filter  can  be  derived  from 
the  point  of  view  of  minimizing  the  variance  of  the  difference 
between  the  predicted  and  actual  noise ,  and  that  this  vari¬ 
ance  must  be  higher  for  any  signal  with  autocorrelation 
different  from  that  of  the  noise. 

III.  Computational  Method 

One  knows  approximately  the  signal  first  arrival  time 
on  all  of  our  seismograms.  In  some  cases  it  1b  directly 
observable,  in  others  one  needs  to  use  travel-time  curves. 

The  autocorrelation  of  the  noise  before  the  first  motion  is 
first  computed.  Prom  this  the  inverse  wavelet  is  computed 
by  the  method  described  in  our  previous  report  Appendix  P 
part  III.  This  is  a  least  iquares  method.  The  length  of 
this  filter  was  chosen  to  be  70  points.  This  is  near  the 
limit  of  computational  feasabllity  of  least  square  proce¬ 
dures  at  the  present  time.  A  method  for  computing  longer 
prediction  operators  was  programmed  but  not  used  because  in 
most  cases  we  did  not  have  a  very  great  amount  of  data 
digitized  before  the  first  motion  and  also  because  exper¬ 
ience  has  shown  that  great  Increases  in  operator  length  do 
not  improve  predictability  proportionately.  Our  data  has 
1/20  second  digitization  intervals,  however,  we  have  dis¬ 
covered  that  our  seismograms  have  little  energy  in  the 
Bpectrum  above  5  ops.  Therefore,  only  alternate  digitized 
points  were  used.  The  resulting  prediction  operator  length 
is  7  seconds. 

Hie  finiteness  of  this  operator  caused  our  actual 
output  to  deviate  from  the  theoretical  output  in  the  follow¬ 
ing  way:  The  operator  cannot  successfully  use  noise  with 
wavelength  of  the  order  of  7  seconds  and  longer  in  predic¬ 
tion,  since  it  is  only  7  seconds  in  length.  Reference  to 
graphs  in  our  previous  reports  indicates  that  5$  to  20# 
of  the  power  in  the  spectrum  may  fall  within  this  range. 
Although  this  low  frequency  is  apparent  in  some  of  the 
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prediction  error  traces ,  the  visual  quality  of  the  records 
is  not  impaired,  however,  due  to  the  very  lowness  of  this 
frequency. 

IV.  Results 

Results  are  presented  in  the  form  of  the  following 
figures.  The  results  are  good  in  every  case  and  sometimes 
quite  remarkable. 
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Figure  VII-5  The  p-wave  is  clearly  located  at  425 
seconds.  This  is  an  example  where  all  but  perhaps  a 
skilled  seismologist  would  not  be  able  to  pick  p  from 
the  original  record,  but  where  it  is  quite  clear  from  the 
prediction  error  record. 
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SECTION  VIII:  Travelling  Auto-Spectra  of  Nuclear  Shot 
Seismograms 

A  travelling  spectrum  is  a  succession  of  spectral 
estimates  of  a  time  function  taken  at  successive  time 
intervals.  Thus  it  is  a  function  of  both  frequency  and 
time.  This  concept,  although  is  is  a  mathematical 
amalgam,  may  be  useful  in  the  analysis  of  non-stationary 
time  series  where  the  successive  spectral  estimates  change 
in  some  physically  meaningful  way. 

It  was  not  oertain  what  could  be  learned  by  taking 
the  travelling  spectra  of  seismograms  of  p-  and  s-waves 
from  nuclear  shots  since  simple  theory  predicts  no  dis¬ 
persion  for  these  phases  in  a  homogeneous  isotropic 
medium.  But  considerable  change  of  waveform  (l.e.  dis¬ 
persion)  is  known  to  occur  in  the  real  earth.  Therefore, 
although  one  has  no  detailed  ideas  of  what  information 
it  might  be  able  to  extract  in  regard  to  nuclear  detec¬ 
tion,  it  was  thought  there  might  be  value  in  computing 
the  travelling  spectra,  especially  since  by  utilizing. a 
special  technique  (Simpson  et  al.,  1961a,  Appendix  J) 
it  was  possible  to  compute  a  travelling  24-point  spectrum 
of  a  typical  Beismogram  on  the  IBM  709  in  the  amount  of 
time  it  takes  to  read  this  sentence.  In  fact,  it  is  too 
easy  to  use  the  computer  to  generate  many  more  numbers 
and  curves  than  are  readily  interpre table.  For  the  first 
investigation  travelling  spectra  was  computed  for  all 
the  digitized  data  which  was  available. 

Since  the  travelling  spectrum  is  a  function  of  two 
parameters,  frequency  and  time,  and  since  our  program  can 
compute  values  almost  as  fast  as  they  can  be  printed, 
there  was  a  significant  problem  in  data  presentation. 

I  took  two  approaches.  The  first  was  to  print  twelve 
numbers  per  line  of  printed  page,  these  being  the  spectral 
amplitude  estimates  scaled  to  a  maximum  of  5,  rounded  to 
an  integer,  and  then  taken  to  the  10th  power.  The  result 


is  intended  to  resemble  13  bar  graphs  running  down  the  page# 
representing  spectral  estimates  at  13  frequencies  as  a 
function  of  time.  Time  of  p-wave  and  s-wave  arrival  is 
indicated.  The  second  approach  to  the  data  presentation 
problem  is  to  make  these  bar  graphs  on  the  scope.  This 
allows  finer  presentation  of  amplitude. 

A  selected  few  of  the  results  are  presented  in  the 
figures*  Some  things  are  notable.  On  Figure  II-3-1  is 
presented  the  travelling  spectra  from  two  nuclear  shots 
over  almost  Identical  paths.  The  spectra  are  similar# 
but  far  from  being  identical.  On  Figure  II-3-5  the  s-wave 
arrival  is  apparent  on  the  travelling  spectra  as  an 
increase  in  high-frequency  energy.  On  Figure  11-3-6  a  phase 
arrival  is  noted  in  which  there  appears  to  be  some  disper¬ 
sion.  This  phase  has  not  yet  been  identified. 


I 


70 


b  LA  j'l  C/L.  0.0.  Km 


.Tt-vur-poi/rrs 

i  ,  s» 


f  %VE 


WAVE 


LOG  A  M 


--Vj/ 


•  r./^LUV,: 

WfCTrttU 


a 

o 

o 

JD 

-Q 

C0 


'tnivelificj  -pcctrci 
f'ji’.Y.c.’L-f  idevrl  icrti 


4rom  +wo  nuckar  sho+s  over 


pa-rKs. 


S 

< 

•#-» 

C9 

O 

Hi 


71 


Figure  VIII-2  Traveling  Spectra  from  Logan  and  Blanca  This  is 
another  type  of  presentation  of  the  information  in  Figure  II-3-1. 
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SECTION  IX  Filtering  for  Signal- to-Noise  Improvement 
Abstract 

A  filter  is  derived  which  can  remarkably  increase 
the  signal-to-nolse  energy  ratio  on  seismic  records.  In 
the  examples  considered  the  ratio  was  increased  by 
factors  of  up  to  about  20.  The  construction  of  this 
filter  is  based  on  assumptions  about  signal  spectra  and 
noise  spectra.  The  filter  distortion,  however,  is  severe 
and  the  method  is  not  expected  to  be  useful  when  applied 
to  first  motion  studies.  Thus  the  method  should  be  use¬ 
ful  for  determining  the  existence  of  very  weak  arrivals. 

A  possible  application  of  this  is  in  the  detection  of 
Lebt's  (Leet  1962)  "lonesome  P"  phase.  This  application 
was  tried  but  results  were  inconclusive  due  to  inadequate 
relevant  digitized  data.  Other  experiments,  perhaps 
less  directly  relevant  to  nuclear  detection,  gave  excellent 
results . 

I .  Introduction 

In  the  previous  subsection  we  have  seen  filtered 
seismograms  in  which  the  signal-to-noise  ratio  was  sub¬ 
stantially  enhanced.  The  filter  in  that  subsection  was 
based  only  upon  a  knowledge  of  the  noise  power  spectrum. 

In  many  geophysical  problems,  some  knowledge  of  the  signal 
may  reasonably  be  assumed.  One  might  make  the  relatively 
weak  assumption  that  the  energy-density  spectrum  of  the 
signal  is  known,  or  one  could  make  the  stronger  assumption 
that  both  amplitude  and  phase  spectrum  (and  thus  the  wave¬ 
form)  were  known.  It  is  advantageous  to  make  the  strong¬ 
est  realistic  assumption  possible  because  -then  the  solution 
filters  are  "tailor-made"  to  the  problem.  It  is  dangerous, 
however,  to  make  strong  assumptions  which  are  not  justi¬ 
fied,  since  we  may  not  know  how  sensitive  our  solutions 
will  be  to  small  deviation  from  the  assumptions.  On  the 
other  hand,  any  sensible  assumption  is  probably  good  if 


the  solution  Is  not  particularly  sensitive  to  deviations 
from  the  assumption. 

II.  Feasibility  Experiment 

In  the  examples  considered  in  this  subsection  we 
assume  knowledge  only  of  the  noise  spectrum  and  the  signal 
spectrum  although  the  method  which  will  be  applied  is 
generally  applicable  to  the  stronger  assumption  of  noise 
spectrum  and  signal  waveform.  The  mathematical  method 
is  to  take  our  assumed  noise  and  signal  spectra  and 
construct  a  filter  which  is  optimum  in  the  Wiener  sense. 

The  details  of  the  method  are  explicitly  developed  in 
Appendices  B  and  C.  The  general  idea  is  that  the  square 
error  will  be  minimized,  error  being  both  l)  filter  out¬ 
put  when  the  only  input  is  noise  and  2)  filter  output 
other  than  signal  when  only  signal  is  input.  It  was 
further  assumed  that  at  a  given  seismic  receiver  noise 
is  present  most  of  the  time  and  signal  is  by  comparison 
rarely  present. 

One  might  wonder  how  sensitive  this  filter  is  to 
small  perturbations  in  the  assumed  signal  and  noise  spectra. 
The  answer  is  that  it  depends  upon  the  spectra.  This  can 
be  seen  by  examining  Figure  HE-'  -1)  in  which  is  displayed 
the  spectra  from  one  of  the  teBt  cases.  The  filter,  as 
might  be  expected,  has  greatest  spectral  components  in 
the  regions  of  high  signal-to-noise  ratio.  It  can  be 
noted  that  high  ratios  at  frequencies  where  both  signal 
and  noise  have  low  energies  do  not  strongly  affect  the 
filter.  Thus  the  filter  seems  to  have  a  sensible  spectrum 
and  although  it  is  peaked  rather  sharply,  it  does  not 
appear  that  any  minor  alteration  in  assumed  signal  and 
noise  spectra  would  cause  major  alterations  in  the  filter 
spectra. 

The  particular  filter  used  in  the  examples  tries  to 
reproduce  the  signal  efter  a  3  second  delay.  To  facilitate 
comparison,  however,  the  time  scale  was  relabled  in  such 


a  way  as  to  remove  the  delay.  Distortion  of  the  signal 
(caused  by  trying  to  suppress  noise)  now  may  cause 
precursers  to  the  signal  as  Boon  as  3  seconds  early.  In 
fact,  the  filter  will  have  considerable  distortion  since 
we  have  set  up  the  problem  bo  that  the  filter  should 
suppress  noise  and  then  we  have  also  said  that  noise  will 
be  the  most  frequent  input.  Thus  the  filter  will  try 
very  hard  to  suppress  noise,  and  much  signal  distortion 
will  almost  always  result.  For  this  reason  the  filter  is 
not  a  good  one  for  first  motion  studies. 

This  part  of  the  experiment  is  based  on  the  follow¬ 
ing  assumptions: 

1)  The  spectrum  of  a  p-wave  signal  from  a  nuclear 
blast  arriving  on  the  LEFT-RIGHT  component  will  be  similar 
to  that  on  the  more  clearly  observable  UP  component. 

2)  The  microseisms  noise  spectrum  does  not  change 
significantly  from  the  minute  before  to  the  minute  after 
a  p-arrival. 

The  results  in  the  particular  examples  studied  which 
are  displayed  in  Figures  (gf2-2  to  3J::  -6)  indicate  that 
these  assumptions  cannot  be  too  bad.  In  them  the  signal 
spectrum  was  determined  from  the  first  25  seconds  of  p-wave 
on  the  vertical  component.  The  noise  spectrum  is  computed 
from  the  horizontal  component  before  the  p-arrival  time. 
(This  time  is  known  from  the  vertical  component.) 

Ill  Detection  Experiment 

In  the  prediction  error  experiment  (Section  VTI  of 
this  report)  one  of  the  prediction  error  filters  increased 
the  signal-to-noise  energy  ratio  to  such  an  extent  that 
the  p-wave  was  easily  recognizable  where  it  had  previously 
required  a  good  deal  of  Imagination  to  recognize  (see 
Figure  k/ti-5) .  Since  this  phase  is  what  Leet  calls 
"lonesome-p"  (more  than  2500  kilometers  distant  and  no 
observable  s -waves  or  surface  waves)  and  its  presence  may 


be  quite  significant  for  nuclear  detection  we  considered 
the  general  problem  of  trying  to  increase  our  ability  to 
detect  just  the  existence  of  a  signal  in  a  very  high 
relative  noise  level.  This  led  to  an  elaborate  mathe¬ 
matical  scheme  written  up  in  SECTION  V.  The  final 
equations  would  be  difficult  to  program  satisfactorily 
using  standard  methods  and  it  was  felt  that  further 
theoretical  Study  would  lead  to  simplifications  both 
theoretically  and  computationally;  therefore,  its  use  is 
not  included. 

The  symmetrical  Wiener-Levlnson  filter, is  quite 
similar  in  concept  and  in  simple  numerical  examples  gave 
similar  numerical  answers.  Furthermore,  one  feels  that 
the  Wiener-Levinson  symmetrical  filter  should  be  able  to 
do  a  better  Job  of  increasing  the  signal-to-noise  ratio 
than  the  prediction  error  filter  because  the  former  is 
derived  from  both  signal  and  noise  information  whereas  the 
latter  is  derived  only  from  noise  information. 

The  essential  assumption  in  this  experiment  is  that 
we  have  some  means  of  getting  knowledge  of  the  lonesome-p 
spectrum.  The  various  possible  means  of  getting  this 
knowledge  represents  a  big  study  in  itself.  In  order  to 
proceed,  we  make  the  following  assumption:  the  spectrum 
will  not  change  radically  from  Logan  to  Blanca  for  similar 
distances  and  similar  paths.  Since  Blanca  was  a  stronger 
blast  than  Logan  it  was  hoped  that  we  would  be  able  to 
find  a  distance  at  which  the  p-phase  could  be  observed  on 
Blanca,  but  ndton  Logan.  Then  we  would  compute  spectrum 
of  the  p-phase  on  Blanca  and  the  spectrum  of  the  noise 
before  Logan  and  construct  a  filter.  This  filter  would 
then  be  applied  to  Logan  in  the  hopes  of  observing  p  on 
Logan.  Unfortunately,  our  available  digitized  data  did 
not  allow  even  this  experiment.  The  closest  approximation 
was  Blanca  3208  km  UP  and  Logan  2111  km  UP.  Unfortunately, 
1)  this  is  nearer  than  the  distances  Leet  specified  for 
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lonesome-p  (2500  km  on  out),  2)  the  distances  may  be 
different  enough  to  cause  a  change  In  the  spectrum,  3)  the 
phase  is  clearly  evident  on  the  Logan  record  even  without 
any  filtering.  The  best  we  can  hope  for  is  that  we  can 
show  improvement  in  the  slgnal-to-nolse  ratio.  Unfortun¬ 
ately,  the  amount  of  noise  digitized  before  the  signal 
arrival  was  so  small  as  the  make  unrealistic  an  estimate 
of  the  improvement  ratio.  Nevertheless,  the  experiment 
was  performed  and  is  depicted  in  Figure  (H-7). 

Better  data  was  clearly  needed. 

IV  Conclusion 

Given  noise  spectra  and  signal  spectra  which  are 
as  different  from  each  other  as  is  typical  with  microseisms 
and  p-waves,  we  can  construct  a  filter  which  substantially 
improves  signal-to-noise  energy  ratio.  Because  of  distor¬ 
tion,  this  filter  is  not  useful  if  a  detailed  study  of  the 
waveform  is  to  be  made. 
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Figure  XC—  2.  E xomple  Of  Wiener  Filtering  Noise  spectre  is  taken 

from  before  the  p-w«ve  on  the  left  component  (top)  and  signal  spectre 

is  taken  from  the  higher  quality  vertical  component  (bottom).  Ihese 
Spectre  ere  used  to  generate  a  filter  which  is  used  on  the  original 

seismogram  %(fop)  to  produce  e  filtered  version  (middle).  Considerable 

Signal  -  to  -  noise  .ratio  improvement  is  noted,  but  the  first 

motion  i*  not  any  more  epparent.  The  exact  time  of  first 

motion  is  known  from  the  vertical  . 
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